######################################
# IMPACT OF COLONIAL RULE ON POLITY SURVIVAL:
#
# Replication Material for:
# Continuity or Change? (In)direct Rule in British and French Colonial Africa
# 
# Carl Mueller-Crepon, 2020
# International Organization
#
# File Description:
# Main analysis of the the survival of precolonial polities during colonial rule
#
# Called from indrule_analysis_all.R 
#
######################################


# Note: 
#' Please note that Stewart's data contains all kinds of polities,
#' not only precolonial ones that have been colonized by either the French or the British. 
#' The data loaded below includes all polities in Stewart's encyclopedia. This is why,
#' in each analysis, the full data is subset to the sample relevant to the analysis:
#' 
# relev.data <- polity.yrs.df[polity.yrs.df$action.type != "decolonization" &
#                               (polity.yrs.df$col.brit == 1 | polity.yrs.df$col.frnc == 1) &
#                               !is.na(polity.yrs.df$start.col) & !is.na(polity.yrs.df$stop.col),]

##########################################
# INIT
##########################################
rm(list = ls())

# Libraries
library(lfe)
library(survival)
library(survminer)
library(stargazer)
library(plyr)
library(ggplot2)
library(raster)
library(rgdal)
library(countrycode)
library(rgeos)

# Paths
fig.path = "figures/polities/"
tab.path = "tables/polities/"

# Output type
output.type <- "latex"

#########################
# FUNCTIONS
#########################

##### LaTeX Functionalities
source("scripts/functions/analysis_functions.R")

# Balance table
source("scripts/functions/balance_table.R")


#################################
# Load and Prepare data
#################################
source("scripts/polities/polities_prepdf.R")



##########################################################
# ESTIMATION SETUP
##########################################################

# ... explanatory variables
ind.vars <- c("col.brit", "year") 
#    polynomial: "x + y + I(x^2) + I(y^2) + I(y*x^2) + I(x*y^2) + I(y^2 * x^2) + I(x * y) "
base.controls <- paste(c("log(popd_1880AD + 1)", "log(dist.coast + 1)", "log(river.dist)", 
                         "log(polity.age+1)"), 
                       collapse = " + ")
ethnic.controls <- paste(c( paste0("v",c(4,5,28,33),".num")), collapse = " + ")
nature.controls <- paste(c("medianaltitude","medianslope","temperaturemean","evapotranspiration",  
                           "precipitation", "ppetratio", "suit", "max.cash.crop" ), collapse = " + ")

# ... fixed effects
strata = " "

# ... SE cluster
se.cluster <- "polity.id"

# ... output specification for stargazer

# ... variables to keep and labels
keep.controls <- c("col.brit",
                   "popd_1880AD", "dist.coast", "river.dist", "polity.age", "year")
keep.labels <- c("British rule", 
                 "Population/km$^2$ (1880, log)", "Distance to coast (log)", "Distance to river (log)", "Polity age (log)", "Year")

# ... notes below tables
notes.p <- "$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01."
latex.notes.add <- function(width = .8){
  paste0("\\parbox[t]{",width,"\\textwidth}{\\textit{Notes:} Cox Proportional Hazard models. Standard errors are clustered on the polity-level. 
         Nature controls consist of median altitude and slope, mean annual temperature, precipitation
and evapotranspiration, the ratio of the two, agricultural suitability, and soils' suitability for cash crop production. 
Ethnic controls consist of the reliance on agriculture and pastoralism, as well as the intensity of agricultural activities. 
         Significance codes: $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01.}")}

latex.notes.add.appx <- function(width = .8){
  paste0("\\parbox[t]{",width,"\\textwidth}{\\textit{Notes:} Cox Proportional Hazard models. Standard errors are clustered on the polity-level. 
         Baseline controls consist of the 1880 population density (logged), the distance to the coast (logged), the age of a polity (loged), and a linear time trend.
 Nature controls consist of median altitude and slope, mean annual temperature, precipitation
         and evapotranspiration, the ratio of the two, agricultural suitability, and soils' suitability for cash crop production. 
         Ethnic controls consist of the reliance on agriculture and pastoralism, as well as the intensity of agricultural activities. 
         Significance codes: $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01.}")}

latex.notes.add.ols <- function(width = .8){
  paste0("\\parbox[t]{",width,"\\textwidth}{\\textit{Notes:} Linear probability models. Standard errors are clustered on the polity-level. 
         Baseline controls consist of the 1880 population density (logged), the distance to the coast (logged), the age of a polity (loged), and a linear time trend.
         Nature controls consist of median altitude and slope, mean annual temperature, precipitation
         and evapotranspiration, the ratio of the two, agricultural suitability, and soils' suitability for cash crop production. 
         Ethnic controls consist of the reliance on agriculture and pastoralism, as well as the intensity of agricultural activities. 
         Significance codes: $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01.}")}

# ... omitted stats
omit.stat <- c("res.dev","ser", "wald", "logrank", "lr")



############################################
# Summary table of all polities
# Appendix Table A5
############################################

source("scripts/polities/polities_table.R")


############################################
# Summary Statistics
# Appendix Table A1
############################################

stub <- "polity.sumstat"

# Select data
sum.df <- polity.yrs.df[polity.yrs.df$action.type != "decolonization" &
                          (polity.yrs.df$col.brit == 1 | polity.yrs.df$col.frnc == 1) &
                          !is.na(polity.yrs.df$start.col) & !is.na(polity.yrs.df$stop.col),]

# Number of polities
length(unique(sum.df$polity.id))
length(unique(sum.df$polity.id[sum.df$countries_cowid == 475]))

# Summary variables
sum.vars <- c("col.brit", "col.frnc", "year",
              "popd_1880AD.ln", "dist.coast.ln", "river.dist.ln", "polity.age.ln",
              paste0("v",c(4,5,28,33),".num"),
              "medianaltitude","medianslope","temperaturemean","evapotranspiration",  
              "precipitation", "ppetratio", "suit", "max.cash.crop" )

sum.labs <- c("British rule", "French rule", "Year",
              "Population (log)", "Distance to coast (log)","Distance to nav. river (log)", "Polity age (log)",
              "Dependence on agriculture", "Dependence on husbandry", "Intensity of agriculture", "Precol. centralization",
              "Altitude (median)", "Slope (median)", "Temperature (mean)", "Evapotranspiration",
              "Precipitation", "Evapotransp. / precipitation", "Suitability for agr.", "Cash crop suitability")

sum.df$popd_1880AD.ln <- log(sum.df$popd_1880AD + 1)
sum.df$dist.coast.ln <- log(sum.df$dist.coast + 1)
sum.df$river.dist.ln <- log(sum.df$river.dist)
sum.df$polity.age.ln <- log(sum.df$polity.age + 1)

# SumStat Table
fileConn<-file(paste0(tab.path, stub, ".tex"))
writeLines(
  stargazer(sum.df[,sum.vars],
            title="Summary of data on lines of succesion",
            notes.align = "l",label=paste0("tab.",stub),align =F,
            digits = 3,  rownames = F,digit.separate = 0,
            covariate.labels = sum.labs,
            type  = output.type), 
  fileConn)
close(fileConn)


############################################
# Balance table
# Appendix Table A15
############################################

# Prepare
west.africa <- c(402:475,483)
coast.west.africa <- c(404,420,433,435,438,452,481,434,437,452,461,471,475)
balance.data <- polity.yrs.df[polity.yrs.df$action.type != "decolonization" &
                                (polity.yrs.df$col.brit == 1 | polity.yrs.df$col.frnc == 1) ,]
balance.data <- balance.data[!is.na(balance.data$polity.id),]
balance.vars <- c("log(popd_1880AD + 1)", "log(dist.coast +1)", "log(river.dist)", "log(polity.age +1)",
                  paste0("v",c(4,5,28,33),".num"),
                  "medianaltitude","medianslope","temperaturemean","evapotranspiration",  
                  "precipitation", "ppetratio", "suit", "max.cash.crop" )
balance.var.labs <- c("Population (log)", "Distance to coast (log)", "Distance to river (log)", "Polity age (log)",
              "Dependence on agriculture", "Dependence on husbandry", "Intensity of agriculture", "Precol. centralization",
              "Altitude (median)", "Slope (median)", "Temperature (mean)", "Evapotranspiration",
              "Precipitation", "Evapotransp. / precipitation", "Suitability for agr.", "Cash crop suitability")
balance.notes <- "$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01. Standard errors are clustered on the polity-level."

# Treatment
treat <- "col.brit "

# Estimate models
model.ls <- list(lapply(balance.vars, function(v){
  felm(as.formula(paste0(v, "~ ",treat," | 0 | 0 | polity.id")),
       data = balance.data)
}),
lapply(balance.vars, function(v){
  felm(as.formula(paste0(v, "~ ",treat," | 0 | 0 | polity.id")),
       data = balance.data[balance.data$countries_cowid %in% west.africa,])
}),
lapply(balance.vars, function(v){
  felm(as.formula(paste0(v, "~ ",treat," | 0 | 0 | polity.id")),
       data = balance.data[balance.data$countries_cowid %in% coast.west.africa,])
}),
lapply(balance.vars, function(v){
  felm(as.formula(paste0(v, "~ ",treat," | fbperp.id | 0 | polity.id")),
       data = balance.data[balance.data$countries_cowid %in% coast.west.africa,])
}))


# Add rows
add.rows <- list(latex.any.row("X-Border FE", c("no", "no", "no","yes")),
                 latex.any.row("Obs", unlist(lapply(model.ls, function(m){m[[1]]$N}))),
                 latex.any.row("British", unlist(lapply(model.ls, function(m){sum(m[[1]]$X[,"col.brit"])}))),
                 latex.any.row("French", unlist(lapply(model.ls, function(m){sum(m[[1]]$X[,"col.brit"] == 0)}))))

# Print
fileConn<-file(paste0(tab.path, "polity.balance", ".tex"))
writeLines(
  balance_table_tex(model.ls = model.ls, treat.var = "col.brit", treat.var.lab = "British",
                    balance.var.label = balance.var.labs,
                    digits = 3, stars = T,
                    label = "polity.balance", col.sep.width = "0pt", title = "Balance test, standardized coefficients", 
                    col.header = c("All", "West Africa", "West African Coast", "X-Border"), 
                    notes = balance.notes, add.rows = add.rows), 
  fileConn)
close(fileConn)


############################################
# Main survival of ever colonized polities
# Table 2, main paper
############################################
stub <- "polity.main"

# ... full specification 
spec <- paste0("Surv(start.col, stop.col, polity.death)", " ~  ", 
               c(paste(c(ind.vars,base.controls), collapse = "+"),
                 paste(c(ind.vars,base.controls,nature.controls), collapse = "+"),
                 paste(c(ind.vars,base.controls,ethnic.controls, nature.controls), collapse = "+")),
               strata,
               "+ cluster(",se.cluster,")")


# Estimation
this.m <- lapply(spec, function(form.str){print(form.str);coxph( as.formula(form.str), 
                                                                 polity.yrs.df[polity.yrs.df$action.type != "decolonization" &
                                                                                 (polity.yrs.df$col.brit == 1 | polity.yrs.df$col.frnc == 1),])})
main.models <- this.m

# Stargazer

# ... Info
add.lines <- list(latex.any.row("Nature controls: ",c("no","yes","yes")),
                  latex.any.row("Ethnic controls: ",c("no","no","yes")))


fileConn<-file(paste0(tab.path, stub, ".tex"))
writeLines(
  stargazer(this.m, se = lapply(this.m, function(m){diag(m$var)^.5}),
            title="British vs. French rule and the demise of precolonial polities: Cox Proportional Hazards",
            dep.var.caption = "End of line of succession",dep.var.labels.include = FALSE,
            keep= keep.controls, 
            covariate.labels = keep.labels,
            notes.align = "l",label=paste0("tab.",stub),align =T,
            add.lines = add.lines, digits = 3, intercept.top = T,intercept.bottom = F,
            omit.stat = omit.stat,
            notes = latex.notes.add(), notes.label = "", notes.append = F,
            type  = output.type), 
  fileConn)
close(fileConn)

# Presentation

# ... Info
add.lines <- list(latex.any.row("Baseline controls: ",c("yes","yes","yes")),
                  latex.any.row("Nature controls: ",c("no","yes","yes")),
                  latex.any.row("Ethnic controls: ",c("no","no","yes")))


fileConn<-file(paste0(tab.path, stub, ".pres.tex"))
writeLines(
  stargazer(this.m, se = lapply(this.m, function(m){diag(m$var)^.5}),
            title="British vs. French rule and the demise of precolonial polities: Cox Proportional Hazards",
            dep.var.caption = "End of line of succession",dep.var.labels.include = FALSE,
            keep= "col.brit", 
            covariate.labels = "British rule",
            notes.align = "l",label=paste0("tab.",stub),align =T,
            add.lines = add.lines, digits = 3, intercept.top = T,intercept.bottom = F,
            omit.stat = omit.stat,
            notes = latex.notes.add(), notes.label = "", notes.append = F,
            type  = output.type), 
  fileConn)
close(fileConn)



###############################
# MAIN SURVIVAL PLOT
# Figure 5, main paper
###############################

# ... main specification 
spec <- paste0("Surv(start.col, stop.col, polity.death)", " ~  ", 
               c(paste(c(ind.vars,base.controls), collapse = "+"),
                 paste(c(ind.vars,base.controls,nature.controls), collapse = "+"),
                 paste(c(ind.vars,base.controls,ethnic.controls, nature.controls), collapse = "+")),
               strata,
               "+ cluster(",se.cluster,")")

# Estimation
this.m <- lapply(spec, function(form.str){coxph( as.formula(form.str), 
                                                 polity.yrs.df[polity.yrs.df$action.type != "decolonization" & 
                                                                 (polity.yrs.df$col.brit == 1 | polity.yrs.df$col.frnc == 1),])})

# Simulate
plot.df <- polity.yrs.df[polity.yrs.df$action.type != "decolonization" & !is.na(polity.yrs.df$start.col) & 
                           polity.yrs.df$start.col == 0 &
                           (polity.yrs.df$col.brit == 1 | polity.yrs.df$col.frnc == 1),]
new.data <- do.call(rbind, lapply(1:0, function(c){
  d <- plot.df
  d$col.brit <- c
  d
}))


# Predictions
m.id <- 1
fit.data <- na.omit(new.data[,all.vars(as.formula(spec[[m.id]]))])
fit <- survfit(this.m[[m.id]], newdata = fit.data)

# Add Survival = 1 at time 0
fit$surv<- rbind(rep(1, nrow(fit.data)),fit$surv)
fit$upper<- rbind(rep(1, nrow(fit.data)),fit$upper)
fit$lower<- rbind(rep(1, nrow(fit.data)),fit$lower)

# Mean survival values
survival.df <- do.call(rbind, lapply(c(1:0), function(c){
  data.frame(Colony = c("French", "British")[c + 1],
             probability = rowMeans(fit$surv[,which(fit.data$col.brit == c)]),
             upper = rowMeans(fit$upper[,which(fit.data$col.brit == c)]),
             lower = rowMeans(fit$lower[,which(fit.data$col.brit == c)]),
             year = c(0:(nrow(fit$surv)-1)),
             stringsAsFactors = F)
}))


#  Survival all obs
survival.all.df <- data.frame(Colony = c("French", "British")[rep(fit.data$col.brit, each = nrow(fit$surv))+1],
                              obs.id = rep(1:ncol(fit$surv), each = nrow(fit$surv)),
                              probability = as.vector(fit$surv),
                              upper = as.vector(fit$upper),
                              lower = as.vector(fit$lower),
                              year = rep(c(0:(nrow(fit$surv)-1)), ncol(fit$surv)),
                              stringsAsFactors = F)

# Compute differences

# ... avergae
survival.diff.df <- survival.df[survival.df$Colony == "British",]
survival.diff.df$probability <- survival.diff.df$probability - survival.df$probability[survival.df$Colony == "French"]
survival.diff.df$Colony <- "Difference"
survival.df <- rbind(survival.df, survival.diff.df)
survival.df$Colony <- factor(survival.df$Colony, levels = unique(survival.df$Colony), ordered = T)

# ... all obs
survival.diff.df <- survival.all.df[survival.all.df$Colony == "British",]
survival.diff.df$probability <- survival.diff.df$probability - survival.all.df$probability[survival.all.df$Colony == "French"]
survival.diff.df$Colony <- "Difference"
survival.all.df <- rbind(survival.all.df, survival.diff.df)
survival.all.df$Colony <- factor(survival.all.df$Colony, levels = unique(survival.all.df$Colony), ordered = T)



#  Plot
png(file.path(fig.path, "polity.surv.main.spagetthi.png"), width = 8, height = 3, unit = "in", res = 600)
g <- ggplot(survival.df, aes(x = year, y = probability, col = Colony)) + 
  
  scale_color_manual(breaks=c("British","French", "Difference"),
                     values = c(`British` = "red", `French`  = "blue", `Difference`  = "black")) +
  geom_step(data = survival.all.df, aes(group = obs.id), alpha = .1, lty = 1, lwd = .5) + 
  geom_step(data = survival.all.df, aes(group = obs.id), alpha = .1, lty = 1, lwd = .5) + 
  xlab("Years since colonization")+ ylab("Survival probability") + ylim(c(0,1)) +
  geom_step(lwd = 1) + geom_step(lwd = 1, lty = 3, col = "white") + #geom_point(pch = 3) +
  theme(legend.position = "none") +
  facet_wrap(~Colony, nrow = 1) + xlim(0,80)
print(g)
dev.off()


############################################
# Main robustness checks
# Appendix Table A10
############################################

stub <- "polity.robust1"

# Main specification
form.str <- paste0("Surv(start.col, stop.col, polity.death)", " ~  ", 
               c(paste(c(ind.vars,base.controls,ethnic.controls, nature.controls), collapse = "+")),
               strata,
               "+ cluster(",se.cluster,")")


# Estimation
this.m <- list()

# ... weighted by colony
this.data <- polity.yrs.df[polity.yrs.df$action.type != "decolonization" &  (polity.yrs.df$col.brit == 1 | polity.yrs.df$col.frnc == 1),]
this.data <- na.omit(this.data[, unlist(lapply(colnames(this.data), function(x){grepl(x,paste(form.str, "countries_cowid"))}))] )
this.data$weights <- gen_weights(this.data$col.brit)
this.m <- c(this.m, list(coxph( as.formula(form.str), data = this.data,
                                weights = this.data$weights)))

# ... weighted by colony
this.data <- polity.yrs.df[polity.yrs.df$action.type != "decolonization" &  (polity.yrs.df$col.brit == 1 | polity.yrs.df$col.frnc == 1),]
this.data <- na.omit(this.data[, unlist(lapply(colnames(this.data), function(x){grepl(x,paste(form.str, "countries_cowid"))}))] )
this.data$weights <- gen_weights(this.data$countries_cowid)
this.m <- c(this.m, list(coxph( as.formula(form.str), data = this.data,
                                weights = this.data$weights)))


# ... stratified by year
this.form.str <- paste0("Surv(start.col, stop.col, polity.death)", " ~  ", 
                                    c(paste(c(ind.vars,base.controls,ethnic.controls, nature.controls), collapse = "+")),
                                    "+ strata(year) ",
                                    "+ cluster(",se.cluster,")")

this.m <- c(this.m, list(coxph( as.formula(this.form.str), data = polity.yrs.df[polity.yrs.df$action.type != "decolonization" &  
                                                                             (polity.yrs.df$col.brit == 1 | polity.yrs.df$col.frnc == 1),])))


# ... additional controls (desease environment)
desease.vars <- c("malaria.tsuit.max", "tsetse.max")
this.form.str <- paste0("Surv(start.col, stop.col, polity.death)", " ~  ", 
                        c(paste(c(ind.vars,base.controls,desease.vars, ethnic.controls, nature.controls), collapse = "+")),
                        strata,
                        "+ cluster(",se.cluster,")")

this.m <- c(this.m, list(coxph( as.formula(this.form.str), data = polity.yrs.df[polity.yrs.df$action.type != "decolonization" &  
                                                                                  (polity.yrs.df$col.brit == 1 | polity.yrs.df$col.frnc == 1),])))


# Stargazer

# ... Info
add.lines <- list(latex.any.row("Robustness check: ",c("empire","colony","stratified","desease")),
                  latex.any.row(" ",c("weights","weights","by year","controls")),
                  latex.any.row("Baseline controls: ",rep("yes",length(this.m))),
                  latex.any.row("Nature controls: ",rep("yes",length(this.m))),
                  latex.any.row("Ethnic controls: ",rep("yes",length(this.m))))


fileConn<-file(paste0(tab.path, stub, ".tex"))
writeLines(
  stargazer(this.m, se = lapply(this.m, function(m){diag(m$var)^.5}),
            title="British vs. French rule and the demise of precolonial polities: Robustness checks",
            dep.var.caption = "End of line of succession",dep.var.labels.include = FALSE,
            keep= keep.controls[1], column.sep.width = "6pt",
            covariate.labels = keep.labels[1],
            notes.align = "l",label=paste0("tab.",stub),align =T,
            add.lines = add.lines, digits = 3, intercept.top = T,intercept.bottom = F,
            omit.stat = omit.stat,
            notes = latex.notes.add.appx(width = .95), notes.label = "", notes.append = F,
            type  = output.type), 
  fileConn)
close(fileConn)


############################################
# Clustering Standard Errors
# Appendix Table A11
############################################

stub <- "polity.clusterse"
cluster.se.spec <- c("countries_cowid","gid")
# ... full specification 
spec <- paste0("Surv(start.col, stop.col, polity.death)", " ~  ", 
               c(paste(c(ind.vars,base.controls,ethnic.controls, nature.controls), collapse = "+"),
                 paste(c(ind.vars,base.controls,ethnic.controls, nature.controls,paste("cluster(",se.cluster,")")), collapse = "+"),
                 paste(c(ind.vars,base.controls,ethnic.controls, nature.controls,paste("cluster(",cluster.se.spec[1],")")), collapse = "+"),
                 paste(c(ind.vars,base.controls,ethnic.controls, nature.controls,paste("cluster(",cluster.se.spec[2],")")), collapse = "+")),
               strata)


# Estimation
this.m <- lapply(spec, function(form.str){print(form.str);coxph( as.formula(form.str), 
                                                                 polity.yrs.df[polity.yrs.df$action.type != "decolonization" &
                                                                                 (polity.yrs.df$col.brit == 1 | polity.yrs.df$col.frnc == 1),])})


# Stargazer

# ... Info
add.lines <- list(latex.any.row("SE clusters: ",c("none","polity","colony","ethnic group")),
                  latex.any.row("Baseline controls: ",c("yes","yes","yes","yes")),
                  latex.any.row("Nature controls: ",c("yes","yes","yes","yes")),
                  latex.any.row("Ethnic controls: ",c("yes","yes","yes","yes")))


fileConn<-file(paste0(tab.path, stub, ".tex"))
writeLines(
  stargazer(this.m, se = lapply(this.m, function(m){diag(m$var)^.5}),
            title="British vs. French rule and the demise of precolonial polities: Standard error clustering",
            dep.var.caption = "End of line of succession",dep.var.labels.include = FALSE,
            keep= keep.controls[1], 
            covariate.labels = keep.labels[1],
            notes.align = "l",label=paste0("tab.",stub),align =T,
            add.lines = add.lines, digits = 3, intercept.top = T,intercept.bottom = F,
            omit.stat = omit.stat,
            notes = latex.notes.add.appx(width = .95), notes.label = "", notes.append = F,
            type  = output.type), 
  fileConn)
close(fileConn)


############################################
# OLS: cross-sectional survival until end of colonial rule
# Appendix Table A12
############################################
stub <- "polity.cross.ols"

# ... prepare
cross.df <- aggregate.data.frame(list(independence = polity.yrs.df$action.type == "decolonization"),
                                 polity.yrs.df[,c("polity","polity.id","countries_cowid","colbrit", "colfrnc")],
                                 FUN = max, na.rm = T)
cross.df <- cross.df[cross.df$colbrit == 1 | cross.df$colfrnc == 1,]
cross.df <- cross.df[cross.df$polity.id %in% unique(polity.yrs.df$polity.id[polity.yrs.df$col.brit == 1 | polity.yrs.df$col.frnc == 1]),]
cross.df <- join(cross.df, polity.yrs.df[polity.yrs.df$start.col == 0 & 
                                           (polity.yrs.df$col.brit == 1 | polity.yrs.df$col.frnc == 1), ],
                 type = "left", by = c("polity.id"), match = "first")


# ... full specification 
spec <- paste0("independence", " ~  ", 
               c(paste(c(ind.vars,base.controls), collapse = "+"),
                 paste(c(ind.vars,base.controls,nature.controls), collapse = "+"),
                 paste(c(ind.vars,base.controls,ethnic.controls, nature.controls), collapse = "+")),
               " | 0 | 0 | ", se.cluster)


# Estimation
this.m <- lapply(spec, function(form.str){felm( as.formula(form.str), cross.df)})


# Stargazer

# ... Info
add.lines <- list(latex.any.row("Baseline controls: ",c("yes","yes","yes")),
                  latex.any.row("Nature controls: ",c("no","yes","yes")),
                  latex.any.row("Ethnic controls: ",c("no","no","yes")))


fileConn<-file(paste0(tab.path, stub, ".tex"))
writeLines(
  stargazer(this.m,
            title="British vs. French rule and the demise of precolonial polities: OLS",
            dep.var.caption = "Reaches independence",dep.var.labels.include = FALSE,
            keep= keep.controls[1], 
            covariate.labels = keep.labels[1],
            notes.align = "l",label=paste0("tab.",stub),align =T,
            add.lines = add.lines, digits = 3, intercept.top = T,intercept.bottom = F,
            omit.stat = omit.stat,
            notes = latex.notes.add.ols(), notes.label = "", notes.append = F,
            type  = output.type), 
  fileConn)
close(fileConn)


############################################
# OLS: Polity Survival year by year
# Appendix Table A13
############################################

stub <- "polity.main.ols"

# ... full specification 
spec <- paste0("polity.death", " ~  ", 
               c(paste(c(ind.vars,base.controls), collapse = "+"),
                 paste(c(ind.vars,base.controls,nature.controls), collapse = "+"),
                 paste(c(ind.vars,base.controls,ethnic.controls, nature.controls), collapse = "+")),
               " | factor(start.col) | 0 | ", se.cluster)


# Estimation
this.m <- lapply(spec, function(form.str){print(form.str);felm( as.formula(form.str), 
                                                                polity.yrs.df[polity.yrs.df$action.type != "decolonization" &
                                                                                (polity.yrs.df$col.brit == 1 | polity.yrs.df$col.frnc == 1),])})


# Stargazer

# ... Info
add.lines <- list(latex.any.row("Year since conquest FE: ",c("yes","yes","yes")),
                  latex.any.row("Baseline controls: ",c("yes","yes","yes")),
                  latex.any.row("Nature controls: ",c("no","yes","yes")),
                  latex.any.row("Ethnic controls: ",c("no","no","yes")))


fileConn<-file(paste0(tab.path, stub, ".tex"))
writeLines(
  stargazer(this.m,
            title="British vs. French rule and the demise of precolonial polities: OLS",
            dep.var.caption = "End of line of succession",dep.var.labels.include = FALSE,
            keep= keep.controls[1], 
            covariate.labels = keep.labels[1],
            notes.align = "l",label=paste0("tab.",stub),align =T,
            add.lines = add.lines, digits = 3, intercept.top = T,intercept.bottom = F,
            omit.stat = omit.stat,
            notes = latex.notes.add.ols(), notes.label = "", notes.append = F,
            type  = output.type), 
  fileConn)
close(fileConn)



##################################
# West Africa: All combined  #####
# Figure 6 main paper
# Appendix Table A16
##################################
stub <- "polity.westafrica.comb"

# Specifications
spec <- paste0("Surv(start.col, stop.col, polity.death)", " ~  ", 
               c(paste(c(ind.vars,base.controls), collapse = "+"),
                 paste(c(ind.vars,base.controls,ethnic.controls, nature.controls), collapse = "+")),
               strata,
               "+ cluster(",se.cluster,")")

# Estimation
this.m <- list()

# ... all West Africa
west.africa <- c(402:475,483)
this.m <- c(this.m, lapply(spec, function(form.str){coxph( as.formula(form.str), 
                                                 polity.yrs.df[polity.yrs.df$action.type != "decolonization" & 
                                                                 polity.yrs.df$countries_cowid %in% west.africa &
                                                                 (polity.yrs.df$col.brit == 1 | polity.yrs.df$col.frnc == 1),])}))




# ... coastal west Africa
coast.west.africa <- c(404,420,433,435,438,452,481,434,437,452,461,471,475)
this.m <- c(this.m, lapply(spec, function(form.str){print(form.str);coxph( as.formula(form.str), 
                                                                           polity.yrs.df[polity.yrs.df$countries_cowid %in% coast.west.africa & 
                                                                                           polity.yrs.df$action.type != "decolonization" &
                                                                                           (polity.yrs.df$col.brit == 1 | polity.yrs.df$col.frnc == 1),])}))



# ... cross border strata
this.strata <- "+ strata(fbperp.id)"
spec <- paste0("Surv(start.col, stop.col, polity.death)", " ~  ", 
               c(paste(c(ind.vars,base.controls), collapse = "+"),
                 paste(c(ind.vars,base.controls,ethnic.controls, nature.controls), collapse = "+")),
               this.strata,
               "+ cluster(",se.cluster,")")
fbperp.west.africa <- c(434,437,452,461,471,475)
this.m <- c(this.m, lapply(spec, function(form.str){coxph( as.formula(form.str), 
                                                 polity.yrs.df[polity.yrs.df$action.type != "decolonization" & 
                                                                 polity.yrs.df$countries_cowid %in% fbperp.west.africa &
                                                                 (polity.yrs.df$col.brit == 1 | polity.yrs.df$col.frnc == 1),])}))



# ... Info
add.lines <- list(latex.any.row("Strata: ",rep(c("--","--", "Border"), each = 2)),
                  latex.any.row("Baseline controls: ",rep(c("yes","yes"), 3)),
                  latex.any.row("Nature controls: ",rep(c("no","yes"), 3)),
                  latex.any.row("Ethnic controls: ",rep(c("no","yes"), 3)))

# ... output
fileConn<-file(paste0(tab.path, stub, ".tex"))
writeLines(
  stargazer(this.m, se = lapply(this.m, function(m){diag(m$var)^.5}),
            title="British vs. French rule and the demise of precolonial polities in West Africa",
            dep.var.caption = "End of line of succession",dep.var.labels.include = FALSE,
            column.labels = collab_w_ruler(column.labels = c("All West Africa", "Coastal West Africa", "x-Border Coastal W. A."), 
                                           column.separate = c(2,2,2), trim = 0),
            column.separate = c(2,2,2),
            keep= keep.controls[1], 
            covariate.labels = keep.labels[1],
            notes.align = "l",label=paste0("tab.",stub),align =T,
            add.lines = add.lines, digits = 2, intercept.top = T,intercept.bottom = F,
            omit.stat = omit.stat,column.sep.width = "-2pt",
            notes = latex.notes.add.appx(.95), notes.label = "", notes.append = F,
            type  = output.type), 
  fileConn)
close(fileConn)

# Plot ###

# ... get coefficients
plot.df <- data.frame(do.call(rbind,lapply(c(list(main.models[[1]]), list(main.models[[3]]), this.m), function(m){
  x <- c(coef = m$coefficients[1], se = diag(m$var)[1]^.5)
  names(x) <- c("coef","se")
  x
})), stringsAsFactors = F)

# ... Add names
plot.df$name <- rep(c("Full sample","All West Africa", 
                      "Coastal West Africa", "x-Border Coastal W. A."), each = 2)
plot.df$name <- factor(plot.df$name, levels = unique(plot.df$name), ordered = T)
plot.df$spec <- rep(c("Baseline", "Nature & ethnic controls"), 4)
# plot.df$spec <- factor(plot.df$spec, levels = unique(plot.df$spec), ordered = T)

# ... sort by name
plot.df$pos <- as.numeric(plot.df$name) + rep(c(-.15, .15), 4)

# ... plot
png(paste0(fig.path, stub,".png"), width = 7, height = 3, unit = "in", res = 600)
g <- ggplot(plot.df, aes(x = pos, y = coef, color = spec)) + geom_point() + 
  geom_errorbar(aes(ymin = coef - 1.96*se, 
                     ymax = coef + 1.96*se), 
                 width=.1) +
  geom_hline(yintercept = 0, color = "darkgrey", lty = 2) +  
  scale_x_continuous(breaks=as.numeric(plot.df$name), labels=plot.df$name) + 
  ylab("Coefficient of 'British rule'") + xlab(NULL) +
  labs(color = "Specification") +
  theme(legend.position = "top")
print(g)
dev.off()


############################################
# Colony - Jackknife
# Appendix Figure A9
############################################

stub <- "polity.coljack"

# ... full specification 
spec <- paste0("Surv(start.col, stop.col, polity.death)", " ~  ", 
               c(paste(c(ind.vars,base.controls), collapse = "+"),
                 paste(c(ind.vars,base.controls,nature.controls), collapse = "+"),
                 paste(c(ind.vars,base.controls,ethnic.controls, nature.controls), collapse = "+")),
               strata,
               "+ cluster(",se.cluster,")")

# ... colonies / countries
colonies <- na.omit(unique(polity.yrs.df$countries_cowid[(polity.yrs.df$col.brit == 1 | polity.yrs.df$col.frnc == 1) & 
                                                           polity.yrs.df$action.type != "decolonization" & polity.yrs.df$ever.colonized]))
colonies <- c(999, colonies)

# ... estimation
this.m <- lapply(spec, function(form.str){
  lapply(colonies, function(c){
    coxph(as.formula(form.str), polity.yrs.df[(polity.yrs.df$col.brit == 1 | polity.yrs.df$col.frnc == 1) & 
                                                polity.yrs.df$action.type != "decolonization" &
                                                polity.yrs.df$countries_cowid != c,])
  })
})

# ... get coefficients
plot.df <- data.frame(do.call(rbind,lapply(unlist(this.m, recursive = F), function(m){
  x <- c(coef = m$coefficients[1], se = diag(m$var)[1]^.5)
  names(x) <- c("coef","se")
  x
})), stringsAsFactors = F)

# ... Add names
plot.df$countries_cowid <- rep(colonies, length(spec))
plot.df$name <- countrycode(as.character(plot.df$countries_cowid), origin = "cown", destination = "country.name.en")
plot.df$name[grepl("Tanzania", plot.df$name)] <- "Tanzania"
plot.df$name[plot.df$countries_cowid == 999] <- " Baseline"
plot.df$spec <- rep(c("Baseline","+ nature controls", "+ nature & ethnic controls"), each = length(colonies))
plot.df$spec <- factor(plot.df$spec, levels = unique(plot.df$spec), ordered = T)

# ... sort by name
plot.df <- plot.df[order(plot.df$name, plot.df$spec),]
plot.df$pos <- rep(c(length(colonies):1), each = length(spec))

# ... plot
png(paste0(fig.path,stub,".png"), width = 7, height = 5, unit = "in", res = 600)
g <- ggplot(plot.df, aes(y = pos, x = coef)) + geom_point() + 
  geom_errorbarh(aes(xmin = coef - 1.96*se, 
                     xmax = coef + 1.96*se), 
                 height=.1) +
  geom_vline(xintercept = 0, color = "darkgrey", lty = 2) +  
  scale_y_continuous(breaks=rev(unique(plot.df$pos)), labels=rev(unique(plot.df$name))) + 
  facet_wrap(~ spec, nrow = 1) + xlab("Coefficient of 'British rule'") + ylab(NULL) 
print(g)
dev.off()





#################################
# Ruler Survival
#################################
source("scripts/polities/analysis_rulers.R")


#################################
# Descriptive plots and maps
#################################
source("scripts/polities/polities_descrplotmaps.R")




